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ABSTRACT 

We report the detection of '^CO 7 = 6^5 emission from the nucleus of the starburst galaxy NGC 253 with 
the redshift (z) and Early Universe Spectrometer (ZEUS), a new submillimeter grating spectrometer This is 
the first extragalactic detection of the '^CO J = 6 —* 5 transition, which traces warm, dense molecular gas. We 
employ a multi-line LVG analysis and find k, 35% -60% of the molecular ISM is both warm (T ^ 1 10 K) and 
dense (hh, ^ 10"* cm""*). We analyze the potential heat sources, and conclude that UV and X-ray photons are 
unlikely to be energetically important. Instead, the molecular gas is most likely heated by an elevated density 
of cosmic rays or by the decay of supersonic turbulence through shocks. If the cosmic rays and turbulence are 
created by stellar feedback within the starburst, then our analysis suggests the starburst may be self-limiting. 

Subject headings: galaxies: individual(NGC 253) — galaxies: ISM — galaxies: nuclei — galaxies: starburst 
— ISM: molecules — submillimeter 



1. INTRODUCTION 

NGC 253 is a nearby (d « 2.5 Mpc; Mauersberger et al. 
1996), highly inclined (i w 78°; Pence 1981) Sc galaxy un- 
dergoing a nuclear starburst (Rieke et al. 1980). The IR lu- 
minosity in the central 30" is 1.5 x 10'" (Telesco & 
Harper 1980), and the mid-IR morphology indicates most of 
this emission arises from a 7" region centered within 1" of the 
2 /im nucleus (Telesco et al. 1993). Near-IR images show a 
prominent stellar bar (Scoville et al. 1985) which likely plays 
a major role in channeling gas into the central starburst (Peng 
et al. 1996; Das et al. 2001). Estimates of the gas mass in the 
centi-al 300 pc range from 0.4-4.2 x 10** M© (Ki-ugel et al. 
1990; Mauersberger et al. 1996; HaiTison et al. 1999). 

The 7 = 6^5 and 7 = 7^6 transitions of CO arise from 
states with energy levels 1 16 K and 155 K above ground, and 
are thus sensitive probes of the warm molecular gas found 
in regions of massive star formation. Observations of these 
lines in the starburst nucleus of NGC 253 (Harris et al. 1991; 
Bradford et al. 2003; Bayet et al. 2004; Gusten et al. 2006) 
have indeed shown that much of the molecular gas is highly 
excited, although a consensus on the details of the excitation 
have yet to be reached (Gusten et al. 2006). Using a multi- 
line excitation analysis, Bradford et al. (2003, hereafter B03) 
find the central 180 pc contain a large mass (2-5 x 10^ M0) 
of warm (T ^ 120 K), dense {n^^ ~ 4.5 x lO'' cm"-') molec- 
ular gas, most likely heated by cosmic rays injected into the 
ISM by the many supernovae (~ 0.1 yr"'; Ulvestad & An- 
tonucci 1997). This model finds large optical depths in the 
mid-J '^CO lines, and thus predicts bright '^'CO emission. To 
further constrain the excitation and energetics of the molecu- 
lar gas we observed the '^CO 7 = 6^5 transition, which pro- 
vides a strong constraint on the '^CO 7 = 6-^5 opacity. This 
is the first extragalactic detection of the '^CO 7 = 6^5 tran- 
sition, and the first detection of any '^CO transition greater 
than 7=3^2 from beyond the Magellanic Clouds. 
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2. OBSERVATIONS 

We observed '^CO 7 = 6^ 5 (433.56 ^m), 7 = 7^6 
(371.65 Aim), i^CO 7 = 6^5 (453.50 ^im), and the [CI] ^Pi 
^Pi fine-sti-ucture line (370.41 ^m) towai'd NGC 253 in De- 
cember 2006 with ZEUS (Stacey et al. 2004) at the Caltech 
Submillimeter Observatory (CSO) on Mauna Kea, Hawaii. 
ZEUS is a direct-detection grating spectrometer providing a 
slit-limited resolving power of A/AA ^ 1000 across the 350 
and 450 /im telluric windows. It currently utilizes a 1 x 
32 semiconductor bolometer array oriented along the disper- 
sion direction, with the pixel size approximately matched to 
a spectral resolution element. A pair of bandpass filters cen- 
tered at 350 /im and 450 /im are mounted directly in front of 
the detector array, such that the system simultaneously pro- 
vides a 16 pixel spectrum in both windows. The bandwidth is 
sufficiently large to simultaneously observe '^CO 7 = 7^6 
and [CI]. 

We obtained absolute spectral calibration with observations 
of Orion (BN-KL), and flux calibration with observations of 
Saturn, which was assumed to have brightness temperatures 
of 116 K, 118 K, and 97 K at 434 /im, 453 //m, and 371 
/tm, respectively (Hildebrand et al. 1985; Orton et al. 2000). 
We observed the four lines over the course of three nights in 
good submillimeter weather, with toisghz = 0.04-0.06. A 
zenith opacity was obtained from both T225 ghz and T350 us- 
ing the CSO atmospheric transmission model, and the mean 
of the two values was used to calculate the transmission to the 
source. Observations of NGC 253 were centered at R.A. = 
00''47"33?2, decl. = -25°17'18" (J2000.0), and small maps 
in the '^CO lines verified that the beam was centered within 
4" of the CO emission peak. We used total power maps of 
Uranus to measure the FWHM of the beam to be 11" at 434 
/tm and 453 fim and 10" at 371 /im. All data were obtained 
by chopping and nodding the telescope with a 30" throw. The 
specti-a of I^CO 7 = 6-^5, '^^CO 7 = 6 -> 5, and the '^CO 
7 = 7^6 and [CI] pair shown in Figure 1 represent total inte- 
grations times of 6, 70, and 5 minutes, respectively. A linear 
baseline is removed from all spectra, and the integrated inten- 
sities are listed in Table 1 . In addition to the nuclear spectra 
we obtained a simultaneous map in the '^CO 7 = 7 ^ 6 and 
[CI] lines, which will be presented elsewhere (T. Nikola et al. 
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Fig. 1.— Top: Spectra of ^^CO 7 = 6 ^ 5 and "CO 7 = 6^5 (scaled 
by X 10). Bottom: Spectmm of '^CO 7 = 7 -+ 6 and [CI] ^ 'Pi with a 
rogue pixel removed near the center. The velocity scale is referenced to '^CO 
/ = 7->-6. 

2008, in preparation). 

3. RESULTS 
3.1. LVG Model 

To examine the CO excitation we assemble the lower-J Une 
intensities from the Uterature, and correct all measurements 
to a common 15" beam. The ^^CO J = 4 —>-3 intensity is ob- 
tained from Glisten et al. (2006), and the 7 = 3 ^ 2 and lower 
ttansitions are taken from Harrison et al. (1999). The *^CO 
/ = 6 ^ 5 map of Bayet et al. (2004) is used to correct the 
intensities measured here to a 15" scale, and when necessary 
the intensities obtained from Harrison et al. (1999) are cor- 
rected using power law interpolations as outiined in BOB. The 
line intensities used in our analysis are listed in Table 1. 

To quantitatively analyze the CO line SED we employ a 
large velocity gradient (LVG) model, in which the excitation 
and opacity of the CO are determined by a gas density (nu^), 
kinetic temperature (Tkin), and CO abundance per velocity 
gradient ([CO/H2\/dv/dr). We use an escape probabihty for- 
malism with /3 = ( 1 - e~^) / r, derived for a spherical cloud un- 
dergoing uniform collapse (Castor 1970; Goldreich & Kwan 
1974). The source is assumed to contain a large number of 
these unresolved clouds, such that the absolute line intensi- 
ties are proportional to a beam-averaged CO column density 
(Nco)- We increase the CO-H2 collisional rate coefficients 
from Flower (2001) by 21% to account for collisions with He 
(B03, and references therein), and fix the H2 ortho/para ratio 
at 3. The CO abundance is set to [CO/H2] = 8.5 x 10"^ (Fr- 
erking et al. 1982) and the isotopologue abundance ratio to 
[i^CO/i^CO] = 40 (Henkel et al. 1993). 



TABLE 1 
Integrated CO Line Intensities 



Transition 


Beam 


I 


I 


a 




["] 


[Kkms"'] 


[ergs s"' cm"" sr 


[%] 


i^CO(6->5) 


11 


854 ± 20 


2.89 X Kr* 


30 




15 


573 


1.94 X 


30 


i2C0(7->6) 


10 


694 ± 30 


3.73 X IQ-* 


30 




15 


41ga,b 


2.25 X lO"'* 


30 


i3CO(6-^5) 


11 


65 ± 5 


1.97 X 10"^ 


30 




15 


43 


1.29 X 10-^ 


30 


i^CO(l^O) 


15 


1105 


1.73 X 10-" 


22 


i2C0(2^1) 


15 


1500 


1.88 X 10-5 


16 


i2C0(3-»2) 


15 


1134 


4.81 X 10-5 


14 


i2C0(4->3) 


15 


1040 


1.04 X 10-^ 


15 


"CO(l^O) 


15 


94 


1.29 X lO-'' 


20 


"C0(2^1) 


15 


113 


1.24 X 10-« 


12 


"CO(3^2) 


15 


117 


4.33 X 10-« 


14 


Note. — Intensities of 




5,7 = 7^6, and 


3C0 7 = 6 



from this work, and lower-I intensities from Gilsten et al. (2006) and Harrison 
et al. (1999). All measurements have been corrected to 15" using a method- 
ology described in § 3.1. Intensities measured here have statistical errors 
listed in column 3, and total uncertainties (30%) dominated by systematic 
uncertainties in the high frequency sky transmission. 

^10% larger than the intensity reported by Gilsten et al. (2006), and well 
within the calibration uncertainties. 

''B03 report an intensity calculated for a source which couples to half the 
power in the Gaussian main beam, rather than to the full main beam as is done 
here. Reducing their value by the corresponding factor of 2 yields an intensity 
10% larger than measured here, well within the calibration uncertainties. 

3.2. High Excitation Component 

As for B03, we find that any single set of LVG model pa- 
rameters capable of producing the mid-J emission underpre- 
dicts the / = 2 — > 1 and 7=1^0 intensities, necessitating 
the adoption of a two component model. Glisten et al. (2006) 
find that a single component can produce the '^CO intensities 
and the 7=3^2 and lower transitions of '^CO, but such a 
model would not account for the bright '^CO 7 = 6 — » 5 enus- 
sion measured here. We begin by using the 7 = 3 — > 2 and 
higher transitions to constrain the high excitation component, 
and then introduce a low excitation component to account for 
the excess 7 = 2^1 and 7 = 1 — > emission. 

We calculate a four-dimensional grid of model CO Une 
SEDs, varying nyi^, Tidn, dv/dr, and A'co over a large vol- 
ume of parameter space. Comparing these model calculations 
to the observed mid-J CO line intensities, we find solutions 
giving Xreduced ^ ^ f^"" valucs of dv/dr > 3 km pc"'. In 
Figure 2 we plot the values of and Tkin giving the best fits 
for velocity gradients in the range dv/dr = 3-320 km s"' pc"^ . 
Over the modeled range of dv/dr these values change by an 
order of magnitude or more, so to further restrict parameter 
space we must apply prior constraints to dv/dr and Tkin. 

In the LVG approximation the velocity gradient is produced 
by large-scale systematic motion, but for a self-gravitating 
cloud in virial equilibrium we can approximate dv/dr w 3.1 

kms"' pc"' ^nHz/lO'* cm'^ (Goldsmith 2001). Allowing that 
dv/dr may be larger due to the presence of an additional stel- 
lar mass density or a high-pressure intercloud medium (B03), 
we set an upper limit ~ 10 times larger at dv/dr < 40 km 
S-' pc"'. Tkin is restricted by the results of Rigopoulou et al. 
(2002), who conclude that the bulk of the warm molecular 
gas ttaced by H2 rotational transitions lies at 7" = 195 K. As 
the mid-J CO transitions arise from lower energy states than 
those producing the H2 rotational lines we expect the mid-J 
CO emission to trace a cooler component, and therefore re- 
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quire T^in < 200 K. With these two upper Hmits the velocity 
gradient is effectively restricted to dv/drK, 7-40 km s"' pc~\ 
with corresponding limits to and (Figure 2). 

To quantify the allowed ranges of the model parameters we 
adopt a Bayesian formalism and calculate a posterior prob- 
ability density function for each parameter (cf. Ward et al. 
2003). We assume a prior expectation of uniform proba- 
bihty per logarithmic interval for each parameter, subject to 
the upper limits on dv/dr and 7id„ imposed above. We find 
103.8- 1041 cm-^ Tkin = 80-200 K, and the thermal 
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pressure is P/ks = 0.8-1 .4 x 10 cm K. The beam-averaged 
CO column density of this warm component is well con- 
strained to be Nco = 1.7-2.2 x 10'^ cm~^, giving an asso- 
ciated H2 mass of = L2- 1.6 x 10^ Mq in the central 
180 pc. We take as our benchmark model the best fit solution 
obtained by fixing dv/dr = 20 km s"' pc"' and plot it over the 
data in Figure 2. 

3.3. Low Excitation Component 

The residual / = 2 ^ 1 and / = 1 ^ intensities from the 
benchmark model can be produced by a broad range of low 
excitation components with 7]dn 40 K and n^^ ^ 1 0^ -10^° 
cm~^, contributing a beam-averaged CO colunm density of 
A^co = 1-5-3.2 X 10'^ cm"^. We therefore estimate the central 
180 pc contain an H2 mass of Mh^ « 2.9 x 10^ M©, 35%- 
60% of which is in a warm {T^^i ~ 1 10 K), dense (whz ~ 10"* 
cm~^) phase. 

3.4. Comparison with Atomic Gas 

Carral et al. (1994) detect 158 ^m [CII] fine-structure fine 
emission and enoission from other ionized and neutral gas 
tracers toward the central 45" of NGC 253. Based on a com- 
bined HII region and photodissociation region (PDR) model 
of the line and far-IR continuum flux, they estimate an atomic 
PDR mass of Mm = 2.4 x lO*" Mq. Scaling from the '^CQ 
y = 1 ^ morphology (PagUone et al. 2004) we estimate half 
of the PDR emission arises from the central 15", givingMei » 
1.2 X 10^ M0 in the central 180 pc, a factor of « 12 smaller 
than the warm molecular gas mass. PDR models generally 
predict comparable amounts of warm molecular and atomic 
components, so we conclude that the bulk of the warm molec- 
ular gas is not UV-heated gas associated with PDRs. In the 
next section we explore alternative mechanisms for heating 
the molecular gas. 

4. DISCUSSION: WHAT HEATS THE GAS? 

4.1. X-Rays 

Due to their smaller cross sections. X-ray photons pene- 
trate more deeply than UV photons into clouds and heat a 
larger volume of the molecular gas. In this section we con- 
sider whether an X-ray Dominated Region (XDR) can pro- 
duce a warm molecular gas mass significantly in excess of the 
warm atomic gas mass. 

Meijerink & Spaans (2005) present the thermal and chem- 
ical structure of four model XDRs, with combinations of low 
or high density («h = 10^ 10^^ cm~^) and low or high inci- 
dent X-ray flux (Fx = 1 -6, 160 ergs s"' cm~^). For each model 
they plot the gas abundances and temperature as a function of 
depth into the cloud, from which we calculate the total col- 
unm densities of warm C"^ and CO. The 158 /zm [CII] transi- 
tion used to trace the warm atomic component arises from a 
state 91 K above ground, so we include only C"^ warmer than 
91 K. We include all CO warmer than 80 K, the minimum 
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Fig. 2. — Results of LVG analysis. Top: Best fit values of n^^ (solid) and 
7iu„ {dashed) as a function of dv/dr, with error bars showing ±1<t range of 
the posterior probability density functions. The dotted lines show the upper 
limits to Tjcu, and dv/dr. Bottom: Integrated line intensities along with the 
benchmark high excitation model (dashed), a representative low excitation 
component (dotted), and the sum (solid). 

temperature allowed by our CO excitation analysis. The large 
observed ratio of NcoINc* w 1.7 (corresponding to M-nJMm 
« 12 as discussed in § 3.4) can only be produced by the high 
density (mh = 10^^ cm"^), high flux (fx = 160 ergs s"' cm"^) 
model. 

The model XDRs use densities which are an order of mag- 
nitude larger or smaller than the value of nn = 2nH2 ~ 10^^ 
cm"-' indicated by our LVG analysis. Interpolating between 
the high density, high flux model and the lower density mod- 
els, we estimate an XDR with hh = 10^'^ cm~^ will match the 
observed Ato/A^c+ ratio only if Fx > 10 ergs s"' cm"^. How- 
ever, such an XDR will produce an [OIJ 63 /im/[CIl] 158 /im 
ratio more than ^ 20 times larger than observed (Carral et al. 
1994; Meijerink et al. 2007), and consequently we rule out an 
XDR as a potential source of the mid-J CO emission. 

4.2. Cosmic Rays 

It is generally accepted that low energy cosmic rays con- 
trol the thermal and chemical balance in the UV-shielded in- 
ner cores of Galactic molecular clouds (Goldsmith & Langer 
1978). B03 estimates that the high supernovae rate in the nu- 
cleus of NGC 253 results in a cosmic -ray ionization rate 750 
times higher than in the Galactic plane, and that these cos- 
mic rays deposit (5-18) x 10"^^ ergs s"' per H2 molecule in 
the molecular gas. By summing the integrated intensities pre- 
dicted by our benchmark LVG model over all rotational tran- 
sitions, we estimate a warm molecular gas mass of « 1 .4 
X 10^ Mq produces a CO luminosity of Zco « 1.6 x 10'' L©, 
corresponding to a specific cooling rate of w 6.9 x 10"^^ ergs 
s"' per H2 molecule. As this cooling rate matches the heat- 
ing rate estimated by B03, we suggest an elevated cosmic-ray 
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density resulting from the starburst may provide the origin of 
the warm molecular gas. 

4.3. Shocks 

The molecular gas in the center of NGC 253 shows ev- 
idence of shock-driven chemistry, including the large gas 
phase abundance of silicon (Carral et al. 1994; Garcfa-Burillo 
et al. 2000), and the general chemical similarity to shock- 
dominated molecular clouds in the Galactic center (Martm 
et al. 2006). In this section we consider whether the mid-J 
CO emission may arise from shock-heated gas. 

The emission from C-shocks is modeled by Draine et al. 
(1983) and Draine & Roberge (1984). In addition to the ro- 
tational transitions of CO, the dominant coolants are the H2 
rovibrational transitions and the 63 /im [01] fine-structure 
line, and we compare observations of these tracers with the 
model predictions in an attempt to constrain the possible 
shock parameters. Engelbracht et al. (1998) suggest that ther- 
mal emission from shock-heated gas in the central 15" pro- 
duces an H2 luminosity of Ln^ = 1.3 x 10^ L0, summed over 
all infrared rovibrational transitions. Carral et al. (1994) mea- 
sure a 63 /zm [OI] luminosity of L[oi] = 8.8 x lO*" Lq in a 42" 
beam, of which we estimate half arises from the central 15". 
With a total CO luminosity of Leo « 1.6 x 10^ we find 
that the L\iJLco ~ 1 and L[oi]/Lco ~ 3 ratios, as well as the 
general shape of the CO line SED, are reproduced for a low 
velocity (Vshock % 8 km s~^) shock incident on nn = lO''- 10^ 
cm~^ gas. 

The molecular gas may be heated by shocks originating in 
the decay of supersonic turbulence, as is the case in the central 
2 pc of the Galactic center (Bradford et al. 2005). Numerical 
simulations of magnetohydrodynamic (MHD) turbulence by 
Mac Low (1999) show that the conversion of dynamical to 
thermal energy in turbulent gas produces a specific luminosity 
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of 

where Vrms and A are the characteristic velocity and length 
scale of the turbulence. The highest resolution molecular 
maps of NGC 253 utilize a ~ 2" beam (~ 24 pc), so structure 
on sub-parsec scales remains unresolved. However, molecu- 
lar gas in the Galactic center shows clumping on scales down 
to ^ 0. 1 pc (Bradford et al. 2005), so we adopt A = 0. 1 pc. For 
a total molecular gas mass of « 2.9 x 10^ Mq the lumi- 
nosity to mass ratio observed in the primary shock coolants 
is « 0.25 L0/M0, comparable to the value obtained from the 
above equation by setting Vmis = Vshock = 8 km s"'. We con- 
clude that in addition to an elevated density of cosmic rays, the 
dissipation of turbulent energy through low velocity shocks 
can produce the warm molecular gas emitting in mid-J CO 
Unes. 

By keeping a large fraction of the molecular gas warm and 
therefore less susceptible to gravitational instability, cosmic 
rays and the decay of turbulence work to halt the starburst. 
As the cosmic rays are produced in supemovae and the tur- 
bulence may be driven by stellar feedback, we suggest the 
starburst in the nucleus of NGC 253 may be self-limiting. 
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